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We investigate features in the single-particle spectral function beyond the Hubbard bands in the 
strongly correlated normal phase of the Bose-Hubbard model. There are two distinct classes of 
additional peaks generated by the bosonic statistics. The first type is thermally activated Hubbard 
“sidebands”, with the same physical origin as the zero-temperature Hubbard bands, but gener¬ 
ated by excitations from thermally activated local occupation number states. The second class are 
two-particle fluctuation resonances driven by the lattice dynamics. In the unity filling Mott insu¬ 
lator, this takes the form of a localized triplon combined with a dispersing holon. Both types of 
resonances also manifest themselves in the structure factor and the interaction modulation spec¬ 
tra obtained from nonequilibrium bosonic dynamical mean-field theory calculations. Our findings 
explain experimental lattice modulation and Bragg spectroscopy results, and they predict a strong 
temperature dependence of the first sideband, thereby opening the door to precise thermometry of 
strongly correlated lattice bosons. 
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I. INTRODUCTION 

Interacting lattice bosons are an active field of re¬ 
search, spurred by continuous experimental advances in 
cold atom systems in ei- Nowadays, not only are the 
experimental parameters such as interactions and opti¬ 
cal lattice depths well controlled [3], but there has also 
been great progress in the field of cold atom spectroscopy 
SHU, giving access to dynamical quantities such as the 
structure factor and the single-particle spectral function 
EHnj. In the strongly interacting Mott regime, there are 
experimental reports on high-energy features in the ab¬ 
sorption spectra obtained by lattice modulation mm 
and Bragg spectroscopy mm- The features appear be¬ 
yond the first Hubbard resonances at higher multiples of 
the local interaction energy, and they are generally inter¬ 
preted as multiple occupation number fluctuations. How¬ 
ever, a detailed understanding of the underlying physics 
is lacking. 

The lowest Hubbard resonances correspond to quasi¬ 
particle and quasihole excitations, whose dispersion can 
be understood already in mean-field and slave-particle 
approaches [I8l - t23] . The full spectral function and struc¬ 
ture factor of the first Hubbard satellites, including fi¬ 
nite lifetime broadening, have so far only been settled 
in one dimension using numerically exact lattice quan¬ 
tum Monte Carlo (QMC) [24] [25] and density matrix 
renormalization group (DMRG) calculations [26]. How¬ 
ever, there are several indications that spectral features 
beyond the Hubbard bands should exist. They have 
been found in the spectral function at zero temperature 
using the variational cluster approach (VCA) j27J [23] 
and strong-coupling calculations [26], and they are also 
reproduced in the current and kinetic energy suscep¬ 
tibilities calculated by DMRG [29]. The recently re¬ 
newed interest in high-dimensional lattice bosons out- 


of-equilibrium, to realize complex effective Hamiltonians 
including gauge fields [5DH52] and spin-orbit interactions 
1331135) by nonequilibrium driving, and the recent real¬ 
time generalization [36] of bosonic dynamical mean-field 
theory (BDMFT) [37H4T] . call for a systematic investiga¬ 
tion of the positions, origins, and temperature behaviors 
of these high-energy fluctuations. 

In this study, we investigate the nature of the res¬ 
onances of the spectral function beyond the Hubbard 
bands at arbitrary integer filling and nonzero temper¬ 
ature. We also propose an interaction modulation spec¬ 
troscopy experiment, and we show using nonequilibrium 
real-time BDMFT [3B, 021 that it probes the two-particle 
scattering susceptibility. Our calculations explain the 
experimental observations from lattice modulation and 
Bragg spectroscopy HMH], while making additional pre¬ 
dictions regarding the temperature dependence, going 
beyond previous zero-temperature calculations in one di¬ 
mension i nr . 

II. MODEL 

We study the strongly correlated limit of cold 
atoms in the first band of a deep optical lattice, by 
means of the Bose-Hubbard model 00] with Hamiltonian 

h = -j + b J b i) + |X h \ h \ b i b i > c 1 ) 

(ij> » 

where h\ ( b { ) creates (annihilates) a boson on site i, U 
is the local interaction due to two-particle s-wave scat¬ 
tering of the neutral bosonic atoms, and J is the nearest 
neighbor lattice hopping integral, which we take as our 
unit of energy. We further limit the study to the normal 
phase without symmetry breaking on the Bethe lattice 
in the limit of infinite dimensions. 
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In this limit the Bose-Hubbard model is described ex¬ 
actly by BDMFT [3714TT] , in direct analogy with fermions 

[35] . The model has a direct bearing on cold atoms in 
three-dimensional cubic lattices, while still being simple 
enough to allow an understanding of the spectral function 
by analytic arguments. In BDMFT the lattice is mapped 
to an interacting impurity coupled to a dynamic external 
bath with hybridization function A(w). For the Bethe 
lattice, the self-consistency directly relates A to the lo¬ 
cal Green’s function of the impurity G(u>) according to 

A (ui) = J 2 G(uj) [35]. 

The impurity model can be solved exactly for nonzero 
temperatures in imaginary time using the continuous¬ 
time quantum Monte Carlo method (CT-QMC) 30j . 
However, Monte Carlo approaches are not able to resolve 
spectral features at high ui due to the required analyt¬ 
ical continuation [46], nor are they suitable for out-of- 
equilibrium calculations due to the dynamical sign prob¬ 
lem m■ Instead, we apply a zeroth- and first-order 
strong-coupling expansion in A, namely the Hubbard- 
I and non-crossing approximation (NCA), respectively 
[551 45], which enables us to directly calculate real-time 
and real-frequency response functions. While the NCA is 
limited to strong interactions, already in this limit there 
are nontrivial additional single-particle excitations in the 
Bose-Hubbard model, as we will show. 

III. METHOD 

The real-time generalization of BDMFT+NCA to 
nonequilibrium situations has been published elsewhere 

[36] . Here we specialize to the equilibrium case, where the 
equations can be Fourier-transformed to real frequency; 
see Appendices |B] and [C] The NCA amounts to mapping 
the local impurity Fock space to pseudoparticles and then 
performing a first-order resummation of hybridization 
events. On the real-frequency axis, the single-particle 
spectral function A is given by the bubble diagram, 

= ~ (2 lyJZ de ( Tr[d< {e)bd> (e+w)5t] 

-Tr [G<( c )&tG>(e-u,)&]), (2) 

in terms of the pseudoparticle propagator G, correspond¬ 
ing to a concomitant fluctuation in occupation number 
states n and n ± 1 on the impurity. The pseudoparti¬ 
cle self-energy is given by the shell diagrams with a 
forward- and backward-propagating hybridization func¬ 
tion, and it takes the real-frequency form 

S > (w) = J 2 f de (f(e)A(e)[bG > (ijj + e)b^] 

J — OO 

+ [l + f{e)]A{e)[tf&>(u>-e)b]), (3) 

where f(u>) = (e^“ — l) -1 is the Bose distribution func¬ 
tion with inverse temperature /3, and the BDMFT lattice 


self-consistency has been used; for details see Appendix 
[C] The first term in Eq. ([3| describes a particle fluc¬ 
tuation on the lattice involving the occupied density of 
states f(u])A(tj), while the second term is a hole fluctua¬ 
tion with the unoccupied density of states [1 +f (uj)]A(tu). 


IV. RESULTS 
A. Analytcial considerations 

1. Hubbard bands 

The main spectral features of the Bose-Hubbard model 
can be understood already at zeroth order in A. This 
amounts to the Hubbard-I approximation (HIA) [35] 
where the lattice self-energy E(w,k) is approximated 
with the zero-hopping (J = 0) self-energy Ehia(+), 
E(w,k) ft Ehia(w). The self-energy can be determined 
analytically, see Appendix [A] and it takes the form 

£hia(z) = [2 Un{z + n) - U 2 n(n - \)]/{z + n + U) , (4) 

at zero temperature, with a single pole at negative fre¬ 
quencies. Reinsertion into the lattice Green’s function 
G(z, k) w [z + p,— ek — Ehia(.s)] , with the noninteract¬ 
ing single-particle dispersion e^, produces two excitation 
branches with dispersion, 

2?k = £k + U (2 n — 1) — 2/j 

± \jei + 2C/(2n + l)e k + G 2 . (5) 

Expanding to third order in J/U the bandwidths W of 
the two branches are given in terms of the noninteracting 
bandwidth W as W = W(n + 1) and IFn. The upper and 
lower branches are centered at e = U(2n— 1 ± l)/2 — fi 
respectively. These two spectral features are the bosonic 
analog of the Hubbard bands in the fermionic Hubbard 
model corresponding to the local particle and hole exci¬ 
tations | n) -A |n±l). We note that the simple Hubbard-I 
approximation reproduces previous mean-field and slave- 
particle results [T5H55] in the normal phase. 

In the fermionic Hubbard model, the effect of dynam¬ 
ical lattice fluctuations is to modify the shape of the 
Hubbard bands. In the Bose-Hubbard model, however, 
the bosonic statistics does not restrict the occupation of 
a single state, and additional fluctuations become rele¬ 
vant. This is directly evident when studying the unity 
filling n = 1 spectral function using BDMFT+NCA; see 
the upper panel in Fig. [T| The upper and lower Hub¬ 
bard bands are located at w ft ±17/2 with approximate 
bandwidths W ft 2 W and W, respectively, with the non¬ 
interacting bandwidth being W = 4 J, in accordance with 
the Hubbard-I approximation. 
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FIG. 1. (Color online) Bose-Hubbard model single-particle 
spectral function A(w) (green solid line) for J = 1, U = 10, 
p = U/2 and the occupied density of states (blue 

dotted line), at high temperature T = 2 (upper panel) and 
low temperature T = 1/4 (lower panel). Inset: pseudoparticle 
Green’s functions — Im[G^ (w)]/7r for the lowest occupation 
number states n = 0- 3atT = 1/4. 


Insertion into Eq. ([3]) gives 


iV 


2tt 


de 


n{m + l)G>_ 1 (w - e)G> +1 (e) 


+ (n + l)mG/) +1 (w - e)G^ l _ 1 (e) . (10) 


Hence, at unity filling the holon and doublon self-energies 
become 


E>M« j 2 g>M, (11) 

E> (w) » 4 J 2 G> (w) + J 2 3Gq (w - 217), (12) 

where in the last relation we have used that the local 
triplon pseudoparticle Gg is long-lived, i.e., (to) = 
l/(w + ir] — 217), as confirmed by numerical calculations; 
see inset in Fig. |T] (and Appendix [P] for the generaliza¬ 
tion to arbitrary filling n). We see that the holon self¬ 
energy Eq only depends self-consistently on the holon 
Green’s function Gj/. In analogy to the lattice self- 
consistency for the Bethe lattice, this results in a semi¬ 
circular holon pseudoparticle spectral function. The dou¬ 
blon self-energy E^ has a similar dependence on Gj 1 but 
with an additional triplon term containing the holon Gq 
with an effective frequency shift of the triplon plus holon 
energy 2 U + U/2 = 517/2. 


B. Bosonic dynamical mean field theory 


2. Additional spectral features 

Interestingly, in addition to the Hubbard bands, there 
are two additional spectral features above the upper Hub¬ 
bard band at w w 317/2 and 517/2. The origin of these 
two spectral features can, in part, be understood by tak¬ 
ing the zero-temperature limit of the NCA diagrams in 
Eqs. ([2]) and ([3]). In this limit, the lesser pseudoparti¬ 
cle Green’s function reduces to the integer-filling n Fock 
state, 

G < (co) ~ — *27r|n)<5(w)(n|, (6) 

and insertion into Eq. © gives the spectral function as 

- i2nA(co) ss ( n + l)G> +1 (w) - nG/l^-w), (7) 

where G> is the nth occupation number state propaga¬ 
tor. At unity filling n = 1 this simple structure of the 
spectral function shows that the spectral weight at nega¬ 
tive frequencies is due to the holon propagator Go, while 
at positive frequencies it is due to the doublon propa¬ 
gator G 2 . Also the pseudoparticle self-energy E 5, in Eq. 
([3]) simplifies by noting that the occupied and unoccupied 
density of states becomes 

-i2nf{u)A(uj) S3 (n + 1 )G> +1 (w), (8) 

-*2tt[1 + f(u))]A(u) S3 -«(?>_!(-«). (9) 


The full BDMFT+NCA calculation at low tempera¬ 
ture confirms these analytic arguments; see the lower 
panel of Fig. [T| The lower Hubbard band is produced by 
the hole excitation |1) —► |0) described by the holon prop¬ 
agator —Go(— ui) with a semicircular shape and band¬ 
width 47. The positive frequency spectrum is given by 
the doublon propagator 2 G 2 (w) with an almost semi¬ 
circular upper Hubbard band (|1) —► |2)) and a lo¬ 
cal triplon excitation combined with a delocalized holon 
(|1) — ► 13) (E> \h}). We also note that the high-temperature 
peak at w ss 317/2 is not present at low temperature; see 
the lower panel of Fig. [T] 

At nonzero temperature, the local occupation on the 
lattice fluctuates around the mean filling, and the simpli¬ 
fying approximation of the occupied states in the pseu¬ 
doparticle propagator breaks down as both the holon and 
doublon will contribute. However, the doublon is easier 
to activate thermally due to its two times larger band¬ 
width compared to the holon. The significant doublon 
contribution is seen in the occupied density of states of 
the upper Hubbard band in the upper panel of Fig. |Tj 
The doublons also contribute to the spectral function, 
which corresponds to a ±1 particle excitation on the 
ground state. Hence, the thermally activated doublons 
undergo the excitations |2) —>• \2 + (±1)), where the hole 
excitation (|2) —> |1)) modifies the lower Hubbard band 
and the particle excitation (|2) —» |3)) creates a long- 
lived triplon (without a holon this time) yielding the new 




























4 



FIG. 2. (Color online) Bose-Hubbard model real-frequency 
single-particle spectral function A{lS) (green solid line) and 
occupied density of states f(io)A(io) (blue dotted line) at three 
boson filling n — 3, J = 1, U = 20, and /i = 5/7/2 at high 
temperature T = 3 (upper panel), and low temperature T = 
1/4 (lower panel). 

spectral feature at w « 3/7/2. The inverted shape of the 
12) —» 13) peak as compared to the thermal occupation 
of the upper Hubbard band can also be understood as 
the doublons in the ground state are being excited to 
the triplon with the excitation energy U — e/, hence the 
low-energy doublons contribute to the high-energy part 
of the peak. 

The additional spectral features described here, 
namely (i) the zero-temperature double-particle excita¬ 
tion | n) — > \n + 2) (g> \(n — 1 )disp) (triplon with dispers¬ 
ing hole at unity filling) and (ii) thermally activated 
side bands \n + 1) —> \n + 2), are fundamental con¬ 
tributions to the Bose-Hubbard spectral function and 
generalize to any integer filling. Unity filling is a spe¬ 
cial case which is nonsymmetric in particle and hole 
excitations. At higher fillings also a double hole exci¬ 
tation is present, | n) —> \n — 2) g) |(n + 1 )disp), and 
thermally activated sidebands proliferate in both posi¬ 
tive and negative frequency as temperature is increased, 

|n± to) —» |n± (to- f-1)), |?i± (to — 1)) (with to < n); see 
Fig- i for the case of n = 3. 


C. Modulation spectroscopy 

To connect the results for the spectral function to 
the particle number conserving system response, we 
now perform interaction modulation spectroscopy and 
compute interaction and density susceptibilities. Using 


U 2U 3 U 



FIG. 3. (Color online) Interaction modulation energy ab¬ 
sorption spectra (solid lines) compared to the local two- 
particle susceptibility Xa+ bt bb (to) (dash-dotted lines) and den¬ 
sity susceptibility XbH^) (dotted lines) for the temperatures 
T = 0.25, 0.5, 1, and 2 (circles, diamonds, triangles, and 
squares, respectively). Inset: real-time evolution of the in¬ 
teraction U(t) and total energy E(t) for driving frequency 
u) = 10 and T = 2. 


nonequilibrium BDMFT, we calculate the energy time- 
evolution E(t) while modulating the interaction sinu¬ 
soidally, /7(f) =Uq + A U sin(wf), with frequency w and 
amplitude AU/Uq = 0.1, during four inverse hoppings 
tmax = 4/J. The energy absorption rate d t E s {t max ) is 
determined from interpolated stroboscopical energy mea¬ 
surements E s (t), see inset of Fig. [3] Sweeping uj in the 
range U/2 < uj < 7C//2 at different initial temperatures 
T produces the spectra of Fig. [3j displaying three res¬ 
onances directly connected to the spectral function fea¬ 
tures. 

The Hubbard-type particle-hole fluctuation (| 1,1) —► 
10,2)) is centered at u « U. The doublon to holon- 
triplon excitation (| 1, 2) —► |0,3)) located at u sa 2 U + 4J 
is strongly thermally activated and features the same in¬ 
verted distribution as the Hubbard sidebands in Figs. 
jT] and [2] The corresponding peak exhibits the addi¬ 
tional shift 4J, i.e., half the bandwidth of the upper 
Hubbard band at low T which is reduced to w « 2 U 
at higher T. The three singlons to two holons and 
one triplon fluctuation (| 1,1,1) —> |0, 0,3)) is located at 
u> ~ 3/7. It exhibits a weak T dependence, and is the 
analog of the triplon plus dispersing holon resonance in 
the spectral function. In linear response, d t E(t) is di¬ 
rectly related to the interaction susceptibility according 
to d t E s (t ) oc (AZ7) 2 Xf,t&tf,;,(w). We compute the impurity 
Xf<t6tf,b(w) using equilibrium real-time BDMFT+NCA, 
and we find quantitative agreement with the absorption 
spectra apart from a slight red shift; see Fig. [3] To con- 
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nect to the structure factor, we also compute the density broken phase is an interesting question for future inves- 
susceptibility x&t&( w ), which has the same structure as tigations. 

Xbtbtbb(w) apart from a 1/2 reduction of the 2 U and 317 
resonances; see Fig. [3] 
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V. DISCUSSION 

Experimentally, the U and 2U resonances have been 
observed in one dimension using Bragg spectroscopy 
ps m and in one and three dimensions using lattice 
modulation spectroscopy rang, where also a weak 317 
resonance has been reported |14j . In contrast to one¬ 
dimensional calculations at zero temperature [43j , we find 
that in higher dimensions and at finite temperatures both 
217 and 317 resonances are present even at commensurate 
filling. We also explain the peak shape of the 217 reso¬ 
nance and its strong thermal activation, almost two or¬ 
ders of magnitude in Fig. [3] confirming indications from 
restricted basis calculations [ 231 . Additional experimen¬ 
tal work in this direction should be worthwhile, especially 
in exploring and using the 217 resonance for Mott state 
thermometry. The 317 triplon excitation is interesting in 
its own right and in connection with Efimov physics [SS3 
152] . Our findings are also relevant in the context of driven 
experiments, e.g., for generating gauge fields [5D1 - {52] and 
spin-orbit interactions [ 331435 ] . where energy absorption 
at the driving frequency needs to be minimized. It would 
also be interesting to study the fine-structure changes in 
the spectral function induced by the higher bands of the 
optical lattice through effective multibody interactions 
[53l ~ f55] and renormalized hopping [56] [57]. 

Theoretically, spectral features beyond the Hubbard 
bands in the normal phase have been reported in one di¬ 
mension at zero temperature using VCA, DMRG, and 
hopping perturbation theory |26H29] without connect¬ 
ing to experiments, while equilibrium temperature ef¬ 
fects have been considered on the slave-particle mean- 
held level [TO] [23]. Within BDMFT in combination with 
CT-QMC, which is exact for the normal phase in infinite 
dimensions, one cannot expect to observe these features 
because of the limitation of analytical continuation ]58| . 


VI. CONCLUSIONS 

In summary, we have analyzed the nature of the reso¬ 
nances of the single-particle spectral function of the Bose- 
Hubbard model in the symmetric Mott phase. We iden¬ 
tified two classes of fundamental excitations beyond the 
Hubbard bands, namely the double-particle excitation 
(triplon plus dispersing holon at unity filling) and ther¬ 
mally activated sidebands generated by excitations from 
thermally occupied number states. These resonances ex¬ 
plain the features in interaction modulation spectra and 
the structure factor, and they are therefore important for 
the interpretation of experiments rara- How these fea¬ 
tures evolve at weaker interactions and in the symmetry- 
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Appendix A: Hubbard-I approximation 

For the Bose-Hubbard model, mean-field calculations 
in combination with the random-phase approximation 
(MF+RPA) and slave-particle approaches already pro¬ 
duce nontrivial spectral functions ram ns ra- In this 
appendix we rederive these results from a DMFT per¬ 
spective and show that in the normal phase the spectral 
functions obtained via these methods are equivalent to 
the Hubbard-I approximation (HIA) [25], where the lat¬ 
tice self-energy is approximated with the self-energy in 
the zero hopping limit J = 0. Without hopping, the 
Bose-Hubbard model separates into a collection of local 
Hamiltonians H = 7 Un(n — 1) — fin. The real-frequency 
spectral function can then readily be obtained from the 
Lehmann expansion of the local single-particle Green’s 
function, 


G l {u) 


1 y- (n\b\m){m\tf\n) , -fj E „ -0E m] 

2 + ir ) + E n - E ™ 

(Al) 


where |?r) is the occupation number state with energy 
E n = \Un(n — 1) — /in. The Green’s function can be 
decomposed into hole and particle excitation branches 
according to 


Gl{u) — G^tuj) + G^\u}) , 


where 


Mp)r \ _ , 1 V" -PE n {n\b\m){m\tf\n) 
G L ^)-+ Z 2^ e w + *77 + E n — E, 

n (h) (i \ _ 1 „- 0 Erl (n\tf\m)(m\b\n) 

L Z io + ir\ — E n + E r 

nm ' 


(A2) 

(A3) 


We immediately see that the negative frequency hole ex¬ 
citations G^ come with negative spectral weight for 
bosons. 
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For integer filling n and zero temperature Eq. (Al) 
simplifies to 


G l (w) = 


1 


ui + irj — Un + fi lo + ir] — Un + U + fi 

(A4) 

This gives the Hubbard-I approximation of the lattice 
self-energy if we use the noninteracting atomic Green’s 
function Gq(z) = l/(z + /i) and the Dyson equation 
E(z) = Gq 1 (z) — G _1 (z). The analytic form for the 


self-energy E is 

S( 2 ) = 


2U n(z + n) — U 2 n(n — 1) 
z + u + U 


(A5) 


Note that the self-energy has a very simple structure with 
a single pole at z = —U — (e.g., for n = 1 and \i = U/2 

it is located at z = —317/2). 

The lattice Green’s function G(w,k) in the Hubbard-I 
approximation is now obtained as 


G(w, k) 


1 

w + irj + fj, — e k - E(w) ’ 


(A6) 


where, as E is momentum-independent, G(w,k) depends 
only on k through the noninteracting dispersion ek- In¬ 
sertion of the analytic form [Eq. (A51] generates two pole 
branches in G(w, k) with dispersions 


2ek = £k + U (2 n — 1) — 2/ii 

± \Jel + 2U{2n + l)e k + U 2 , (A7) 

corresponding to the nil particle and hole excitations, 
i.e., the upper and lower Hubbard bands, with band cen¬ 
ters at e = U(2n — 1) ± l)/2 — /u,. The bandwidths W 
of these bands as a function of noninteracting bandwidth 
W of e k G {—VF/2, IT/2}, W = e£\w/2)±e£\-W/2), 
is to second order mW/U given by W = W(n + 1) and 
W = Wn for the upper and lower Hubbard band, respec¬ 
tively. Hence the widths scale with the integer filling n 
and exhibit the largest asymmetry at unity filling n = 1. 
The Hubbard-I approximation directly reproduces the re¬ 
sults for the normal phase obtained by MFiRPA and 
from slave-particle theory m DU HD m], and it also 
generalizes trivially to nonzero temperatures. 


Appendix B: BDMFT+NCA equilibrium real-time 
propagation 


numerically easier to obtain high-quality real-frequency 
results by performing a real-time evolution and then 
Fourier transforming the results to real frequency. Us¬ 
ing this scheme, one avoids inverting the Dyson equa¬ 
tion on the real-frequency axis, which requires a care¬ 
ful discretization and handling of Lorentzian broadening 
factors. (In the low-temperature limit this inversion is 
increasingly difficult as the pseudoparticle ground state 
approaches a delta function at zero frequency in the Bose- 
Hubbard model.) 

For the time evolution in equilibrium, only the Dyson 
equations for the retarded component G R (t) and the 
right-mixing component G n (t,r) of the pseudoparticle 
Green’s function G are required, assuming that the Mat- 
subara imaginary-time Green’s function G M (r) is known. 
The Dyson equation for G R (t) takes the Volterra form 


(id t — H)G R (t) - f dt£ R (t-t)G R (t) = 0, (Bl) 
Jo 

with the initial boundary condition G R ( 0 ) = — il and 
the pseudoparticle self-energy £ R . The Dyson equation 
for G n (t,r) also depends on £ R according to 

(id t ~ H)Gr (t,r) — / di£ R (t— t)G^(t,T) = , 

Jo 

(B2) 

where the additional right-hand side Q (t, r) is given by 
the Volterra-type convolution of the right-mixing pseu¬ 
doparticle self-energy E (f, r) and the imaginary-time 
pseudoparticle Green’s function G M (r): 


Q n (t,r)=/ df E n (t,f)G M (r). (B3) 

Jo 

The lesser Green’s function G < (t) is directly given by 
the right-mixing Green’s function at r = 0, G < (t) = 
G n (f, 0). Furthermore, the projection onto the physi¬ 
cal space [35] also gives the (pseudoparticle-specific) re¬ 
lation for the greater component G R (t) = 6(t)G > (t). 
Both lesser and greater components G^(t) are readily 
extended to all times using the antihermicity relation 
= -[G^(f)] t g2J. Hence solving Eqs. {Bl} and 
(B2) determines all Keldysh components of G. 


Static local observables are obtained as direct traces 
over the lesser pseudoparticle Green’s function, 


(O) = *Tr[OG<(0)], 


(B4) 


To calculate real-frequency properties in 
BDMFT+NCA we adapt the real-time out-of- 
equilibrium formulation of Ref. [36] to equilibrium, 
where all response functions are time translation in¬ 
variant, G(t,t') = G(t — if). Using the notation of Ref. 
[42] . this yields a simplified set of equations for the 
pseudoparticle real-time propagation. 

While it is indeed possible to transform the resulting 
equations directly to real frequency it turns out to be 


at relative time t = 0, trivially giving the same re¬ 
sult as the imaginary-time equilibrium calculation ( O) = 
-Tr [OG m (/?)] = fTr[OG n (0,0)] = *Tr[OG<(0)]. 

In NCA, the lesser and greater single-particle Green’s 
functions G^ are obtained as the pseudoparticle bubble 

G^(t) = m[G^{-t)bG*{t)tf}. (B5) 

Note that we only consider the normal phase here with¬ 
out symmetry breaking, ( b) = 0, hence the full Green’s 
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function is equal to the connected Green’s function. The 
form for the single-particle Green’s function can be gen¬ 
eralized to arbitrary susceptibilities \\ is ft) between a 

pair of operators A and B yielding 


XA,B® = m( r ^[G < (-t)AG>(t)B] 

- Tr[G> (~t)B G< (f)i]) , (B6) 


i.e., Xb li ft) corresponds to the retarded single-particle 
Green’s function, \b,bsft) = G R (t). Note that the equal 
operator susceptibility is denoted as Xa = Xa A- 

The retarded pseudoparticle self-energy is obtained 
from the greater component E^(l) = 0(1)S > (1) which 
in turn is given by the two shell diagrams, 


E>(i) = i(A>(t)[b f G > (t)b] 

+ A<(-t)[bG>(t)bft), (B7) 


and the right-mixing component is obtained analogously, 


E n (t,r) = i^Aftt,T)[b f G n (t,r)b\ 

+ A r (T,t)[bG~'(t,T)bft'j . (B8) 


Use of the Bethe lattice self-consistency relation A (t) = 
J 2 G(t), then yields the closed set of BDMFT+NCA 
equations in combination with Eqs. (Bl), (B2), (B3|, 
p5|, dB7l), and (|B8). 


The equilibrium real-time Dyson equations are solved 
using an equidistant time-discretized second-order prop¬ 
agation method. While the convolutions scale quadrat- 
ically with the time step G{N.f) as in the out-of- 
equilibrium method, the time-translation invariance re¬ 
duces the memory scaling from quadratic to linear 0(N t ), 
enabling much longer time calculations, which yield a 
very fine real-frequency resolution after Fourier transfor¬ 
mation. 


Appendix C: Real frequency 


While we perform the numerical calculations in real 
time it turns out that formulating NCA in real frequency 
is a fruitful venue for understanding the physics, es¬ 
pecially the low temperature limit. Hence we Fourier- 
transform the equilibrium real-time equations using 

/ OO 1 /»00 

dt e iuit G(t ), G(t) = — duse~ iuit G{us). 
-oo 2 tt J_ oa 


The single-particle spectral function is obtained from Eq. 
pk as, 


A(u) = --Im[G fl (w)] = ^-[G>(w) - G<(w) 

7 r Z7r 


— r 

( 2^) 2 y _< 


de f Tr[G< (e)b G> (e + oj)bft 

O ' 

- Tr[G < (e)& t G > (e - w)&]) , (Cl) 


and the greater pseudoparticle self-energy E > in Eq. (B7) 
takes the form 


E > (w) = — [ de (A < (e)[bG > (co + e)b^] 

J — OO ' 

+ A>(e)[6tG>(w-e)6]). (C2) 

Now, using the Bethe lattice self-consistency A (to) = 
J 2 G{us) in combination with the relation of the greater 
and lesser components to the single-particle spectral 
function, G > (w) = —i2irA(us)f(u}) and G < (w) = 

— i2'kA{u})[1 + /(<*;)], where f(u) = (e@ u — 1) _1 is the 
Bose distribution function, this can be rewritten as 


E>(w) = J 2 [ de (f(e)A(e)[bG > {u + e)bft 

J — OO ^ 

+ [l + /(e)]A(e)[6 tG>( w -e)6]). (C3) 


These results, Eqs. (Cl I and (C3), correspond to Eqs. 
([2]) and (|3j). 


Appendix D: Analytical results 

In the low-temperature limit of the Mott insulator with 
integer filling n , the NCA relations can be simplified by 
observing that the lesser pseudoparticle Green’s function 
G < , corresponding to the occupied pseudoparticle den¬ 
sity of states, approaches a delta function in real frequen¬ 
cies 


G < (w) ~ —i2Tr\ri)S(ui){n\ . (Dl) 

This follows from Ref. EH where it is shown that the 
energies of the projected pseudoparticles are strictly pos¬ 
itive and the lesser Green’s function is given by G < (w) = 
—i2/G(w)Im[G > (w)], where /g(w) = e -/3aj is the classi¬ 
cal Gibbs distribution function. This suppresses all but 
the local atomic ground state in the zero-temperature 
limit. Insertion into the bubble diagram for the spectral 
function [Eq. ( |C1| )] then reduces this expression to 

- i2nA{ui) « (n + l)G> +1 (w) - nG^^-w). (D2) 

So, in the low-temperature limit, the single-particle ex¬ 
citations are the n — 1 and n + 1 pseudoparticles with 












local interaction energies E n ± \ — E n = U/2. Also the oc¬ 
cupied and unoccupied single-particle states can be ap¬ 
proximated according to 


-i2irf(uj)A(u) « nG^_ 1 (—w), (D3) 

+ f(w)]A(u) « (n + 1 )G> +1 (w). (D4) 


Insertion into Eq. (C3) gives the simplified greater pseu¬ 
doparticle self-energy 




de 


n(m + (w - e)G> + 1 (e) 


-(-(n-(-1 )toG^ +1 (w - e)G^_ 1 (e) , (D5) 


which corresponds to Eq. |To| . In particular the nil 
pseudoparticles obey the relations 


iJ 2 


de 




i (n i l)(n - 1 ) G> +1 (w - 


_ e )Gn( e ) 

e )^n- 2 ( e ) 


(D 6 ) 


and 


^+iM 


u 2 r°° 

27T ^ 


de 


n(n i 2)G>_ 1 (w - e)G> + 2 (e) 


i (n i l ) 2 G> +1 (w - e)G^(e) . (D7) 

This can be si mpl ified using the low-temperature approx¬ 
imation [Eq. (D 1 )] once more, using G>(w) ~ — i2-kS(uj). 
However, we go one step further and also apply the nu¬ 
merically motivated approximation of sharp n ± 2 pseu¬ 
doparticle resonances, i.e., G> ±2 (w) « — i2Trd(aj— [E n ± 2 — 
E n ]), where E n ± 2 — E n = 2 U; see also inset in Fig. 1. 
The local interaction energy E m of the occupation num¬ 
ber state | m) is given by E m = Um{m — l)/2 — /im with 
fixed /i = (2 n — l)U/2. Both of these approximations 
inserted into Eq. (D7) give 


« J 2 n 2 G>-tM + J 2 {n 2 - l)G> +1 (w - 2 U ), 

E> +1 ( W )« J 2 (n + 1) 2 G> +1 M 

+ J 2 n(n + 2)G>_ 1 {u-2U). (D 8 ) 


The first term in these relations, the reappearance of 
the pseudoparticle Green’s function G > ±1 in its own self¬ 
energy £> ±1 , is the dominant self-energy effect. It acts 
similar to the Bethe lattice self-consistency relation, pro¬ 
ducing semicircular spectral functions for the nil pseu¬ 
doparticles centered around w = E n ±i — E n = U/2. 
Their respective bandwidths are directly obtained from 
the prefactors as W n _\ = inJ and W n+1 = 4(n+ 1) J, in 
agreement with the Hubbard-I approximation. The sec¬ 
ond term in the self-energy expressions is shifted up in en¬ 
ergy by 2/7, hence it is located at w = U/2 + 2/7 = 5/7/2, 
and it gives a small correction to the spectral functions. 

Note that the relations for the n — 1 self-energy only 
hold if n > 2. In the special case of unity filling n = 1 
there is no n — 2 Fock state available, and the holon self¬ 
energy simplifies to 


£>(u,)« J 2 G 5 », (D9) 

while the doublon self-energy takes the form 

£> M » 4 J 2 G> (w) + J 2 3G> (w - 2/7). (DIO) 


These are the simplified pseudoparticle self-energy rela¬ 
tions in Eqs. and ( p~2] > . As the holon self-energy 
So only depends on the holon propagator Gq , it fol¬ 
lows immediately that Gq is semicircular and centered at 
u = U/2. By Eq. (D2 1 , — Gg”(— u) directly gives the lower 
Hubbard band in the spectral function A{iS). The dou¬ 
blon self energy T, 2 has the same type of semicircular gen¬ 


erating term 4 J 2 G 2 but also a high-energy triplon plus 
holon correction emerging in terms of the holon propa¬ 
gator J 2 3Gq(lo — 2/7). Hence in the positive frequency 
spectral function, the upper Hubbard band is generated 
by the doublon-doublon self-consistency and the higher 
triplon excitation by the holon correction term; see Fig. 
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